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_  ,  Abstract  *"  '***'  ^ 

%  / 

In  qw  previous  work  {Ik  a  procedure  was  de- 
rithav  fpr  rational  function 

to 


1.2  Differential  Equations  and 
Their  Asymptotic  Solutions 


veloped,  based  on  algor 


r-*—  - r — »  - - - - -  - 

frequency  extrapolation  and  V-norm  averaging 

circumvent  both  the  effects  or  noise,  and  erri 


errors 

due  to  refracted  curved  paths,  for  solving  the  in¬ 
verse  scattering  ^problem  in  ultrasonic  imaging. 


Several  derivations  exist  in  the  literature 
of  partial  differential  equations  which  govern  the 
way  sound  travels  through  B.  For  example,  Che 
Helmholtz  equation  model  is  given  by 


In  this  paper,  we  reverse  the  order  of  the  first 


two  algorithms  with  fespect  to  ID;  that  is,  .we 
first  apply  the  ^/  averaging  procedure  and  then 
carry  out  the  rational  function  procedure,  extrap¬ 
olating  to  infinite  frequency.  An  adaptive  method 
for  choosing  the  best  .  r^t ional  function  expansion 
is  also  employed.  describe  the  underlying 
ideas  of  the  procedure,  and  we  illustrate  examples 
of  reconstruct  ions  in  the  presence  and  absence  of 
Gaussian  noise.  These  results  are  then  compared 
with  results  of  previous  experiments.  The  result¬ 
ing  procedure  works  with  a  much  smaller  signal -to- 
noise  ratio.  An  extension  is  suggested  to  image 
speed  of  sound  and  density. 
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The  Ricatti  model, 
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is  obtained  from  (1.1)  by  applying  the  transform 
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The  Rytov  model  is  obtained  from  (1.3)  by  dropping 
the  term  containing  Vw  • 


1.  Mathematical  Basis  for  the  Algorithms 


A  differential  equation  governing  the  flow  of 
sound  through  B  may  depend  on  other  properties  of 
B,  such  as  the  density  p(r)  and  the  compressibil¬ 
ity  *(r) .  In  general,  such  an  equation  may  take 
the  form 


1.1  The  Physical  Problem 


Au  ♦  F(r,  p(r),  ic(r),  u,  k)  ■  0 


(1.5) 


Let  us  assume  that  sound  at  high  frequency 
penetrates  a  body  B  iMersed  in  fluid.  We  denote 
the  spatial  sound  pressure  and  the  corresponding 
sound  source  by  u  and  ul,  respectively.  The  ve¬ 
locity  of  sound  in  the  fluid  is  given  by  c_ ,  and 
the  velocity  of  sound  in  the  body  by  c(r/.  A 
particular  path  along  which  sound  travels  from  a 
source  point  r  to  a  detector  point  r.  is  denoted 
by  L  if  it  is  a  straight  line  path  and  by  P  if  it 
is  a  curved  path. 


where  A  is  a  differential  operator. 


The  success  of  the  asymptotic  procedure  de¬ 
pends  on  our  knowing  an  asymptotic  expression  for 
some  function  +(u,  k)  which  we  can  compute  or 
measure  for  a  range  of  values  of  k.  For  example, 
+(u,  k)  may  have  the  asymptotic  expression 


♦(u.  k)  -  p(k)[v(«,  p)  ♦  e(k,  «.  p)]  (1.6) 


g(g  FILE  COPY 


.•v 


I 


where  «(k)  is  «  known  function. 


(1.13) 


e(k,  x,  p)  -  o(k  °)  (1.7) 

•s  k  ♦  •,  and  where  o  is  a  positive  number.  In 
addition,  it  is  desirable  to  know  th«t  4(u,  k)  is 
an  analytic  function  of  k  in  some  sector 
{keC:  |arg(k  -  k  ) |  <  d } .  We  can  then  compute 
v(x,  p)  via  the  fhiele  algorithm,  even  though 
c(k,  r,  p)  may  not  be  small  over  the  range  of 
values  of  k  on  which  we  know  4* 


l in  R(k)  •  p! 

k—  2m 

i.e.,  we  have  the  approximation, 


vU.p)  -  P®B.  (1.1*) 

For  example,  for  the  case  of  the  Helmholtz 
model  (1.1),  we  have 


Now,  if  4  is  known  and  can  be  evaluated  for  k 
on  a  finite  ^nte^val,  such  as 

I  •  (k:k  •  -k-k  },  then  we  can  apply  the 
1  min  max 

Thiele  algorithm  to  tne  2m  ♦  l  values 


4(u,k) 


log  u^(r) 


/  /l  ♦  f(r)d»  ♦  0(k*°), 
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♦(u.kj) 
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j  »  0,  1,  ....  2n>; 

k.tl,  k.  distinct. 
J  J 


(1.8) 


to  accurately  estimate  v(«,p).  The  .Thiele  algo¬ 
rithm  (2,  3,  4]  involves  using  p-^  defined  in 

(1.8),  and  then  computing  the  numbers 


(1.15) 

in  which  v  •  J  A  ♦  T  ds,  for  the  esse  of  the 
Rytov  approximation  we  have 


4(u,k) 
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/  f (r)ds  +  0(k~°),  k  ♦  - 
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(1.16) 


k.  -  k. 
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1  pj  +  1  -  PJ 
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0,  1,  ....  2m  -  1 


k.  .  -  k. 

pJ  .  pJ^  ♦  Hi - -L. 

Pi  Pi-2  j  +  1  j 
or  -  or 


"i-l 


"i-1 


j  ■  0,  . .  2m  -  i 

i  «  2,  3,  ....  2m. 


(1.9) 


Some  of  these  numbers  can  be  used  in  s  truncated 
continued  fraction  expansion  of  s  rational  func¬ 
tion  which  interpolates  the  data  p^.  That  is, 
this  rational  function  R(k)  is  given  by  the  con¬ 
tinued  fraction 


R(k)  •  p° 


k-k 


k-k 


2m 


o  o 
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o  o 
p2m  *  p2m-2 


(1.10) 


and  satisfies 


where  in  (1.16),  L  is  a  straight  line  path  join¬ 
ing  r  to  x  .  and  v  *  /  f  ds .  The  examples  of  this 
paper  are  Dated  upon'  using  (1.14)  and  (1.16), 
Eight  significant  figure  accuracy  is  attainable, 
using  ll(m  *  5)  or  13(m  ■  6)  values  for  a  4  mm 
(FWHM)  object  with  no  of  noise.  However,  the 
presence  of  noise  diminishes  this  accuracy  consid¬ 
erably,  although  it  is  possible  to  admit  several 
powers  of  ten  times  as  much  noise  vis  the  use 
of  averaging. 

If  the  influence  of  noisy  data  is  not  unduly 
destructive,  we  can  use  the  above  procedure  to 
recover  two  functions  f(x,p)  and  g(x,p)  provided 
that  a(k)  and  6  >  0  is  known,  and  o  >  0  in  the  ex¬ 
pression  , 

4(u,k)  *  a(k)  |f  ( x,  p)  ♦  k  5[g(x,p)  ♦  0(k  °)]} 

(1.17) 

In  this  case,  we  can  recover  f(x,p)  by  sampling 


0.  1. 


2m 


(1 .11) 


(1.18) 


Once  f(x,p)  has  been  determined  vis  (1.14),  we  can 
It  is  easy  to  see  that  the  rational  function  apply  the  Thiele  algorithm  to 
(1.10)  also  has  the  form 


R(k) 


pL>" 


♦  C,k 


♦  c 


♦  dl  k*=rT 


(1.12) 


♦  d 


where  pP  is  the  last  entry  computed  in  the  table 
(1.9),  end  where  cj  and  d}  (i  •  1,  2,  ...,  m)  are 
the  earns  constants.  Hence,  in  addition  to  (1.11), 
R(k)  also  satisfies  the  relation. 


v  “  k4I^7kT>  '  f(,c-p)l  (!•»*) 

to  get  g(x,p). 


2.  Averaging 

In  [1),  averaging  was  carried  out  after 


r 


application  of  the  Thiele  algorithm.  We  have 
since  found  that  we  can  tolerate  considerably  more 
noise  in  the  data  by  applying  £  averaging  before 
application  of  the  Thiele  algorithm. 

Since  the  data  are  complex  numbers,  actual 
£  averaging  is  too  difficult  numerically.  In¬ 
stead  ,  we  carry  out  an  £.  average  on  both  the  real 
and  imaginary  parts  of  uie  dsta  separately.  The 
process  is  very  simple,  and  works  as  follows. 

Given  a  set  of  rejl  numbers  {a.}*^»  we  want 
to  find  a  real  number  c  which  minimises  the  sum 

S(c)  -  I  |*.  -  c |  (2.1) 

i-1 

The  solution  is 

c*  -  a.  (2.2) 

] 

where  Sj  is  one  of  the  values  of  N°te 

that  for  even  M,  j  may  have  two  solution  values. 
Por  odd  M.  c  is  the  median,  hence,  the  problem  of 
finding  c*  reduces  to  finding  the  minimus  of  the 
values 

»(•;)  -  !  h  -  «jl 

J  i-1  J 

with  respect  to  j. 

3.  An  Adaptive  Feature 

The  method  of  £  averaging  provides  a 
convenient  procedure  for  adaptively  choosing  the 
number  of  terms  in  the  rational  function  .R(k), 
Por  noisy  data  or  for  slowly  changing  p)  with 
respect  to  k;,  a  better  R(k)  is  obtained  by  a 
rational  function  with  a  small  number  of  terms. 

In  our  test  examples  which  follow,  the  method 
of  £  averaging  was  also  applied  to  the  odd  rows 
(i  -*0,  2,  4,  of  the  entries  in  the.p  table. 

For  example,  if  all  of  the  values  of  pJ  in  (1.8) 
are  the  same,  there  is  no  point  in  continuing, 
since  this  value  is  then  also  the  limit  u  k  ♦  * 
of  R(k).  Due  to  round-off  error,  the  criterion 
actually  used  was  to  construct  the  entries  of 
a  particular  row  M2iM#  find  the  £^  average,  c, 
with  respect  to  j  as  described  above,  and  then 
terminate  the  computations  .whenever  the  number 
of  p^  (i  fixed)  such  that  | pJ  -  c |  -6  was  greater 
than1  the  number  of  p?  such  1  that  |  p?  -  e|  >  6  for 
some  small  6.  1  1 

4.  Numerical  Experiments 


Pig,  1.  Digitised  image  of  test  objects  on  a  41  x 
41  grid  with  each  pixel  having  dimensions 
of  .3125  mm  x  .3125  mm.  The  FWHM  for 
each  of  the  Gaussian  distribut ions  is  1.0 
mm  and  their  centers  are  located  at  y  - 
3.75,  0.0,  -2.5,  and  -4.375  «n.  The 

bright  line  through  the  centers  of  these 
objects  indicates  the  location  of  the 
dats  points  for  the  intensity  profile 
which  is  plotted  on  the  right  side. 

where  b  was  chosen  so  that  the  full  width  at  half 
maximum  of  each  object  was  1  m  (which  is  1.67 
wavelengths  at  the  central  frequency  of  2.5  MHz 
with  Che  speed  of  sound  c0  *  1500  m/a)  and  the 
points  r-  were  chosen  at  x  •  0  and  y  -  3.75,  0.0, 
-2.5,  acA  -4.375  mm,  respectively.  All  of  the 
digitised  images  which  we  shall  present  for  this 
test  case  will  use  a  41  x  41  grid  with  each  pixel 
having  dimensions  .3125  mm  x  .3125  mm  and  with  a 
profile  of  the  brightened  central  line  plotted  to 
the  right  of  the  image. 

The  Rytov  model  was  used  to  generate  the 
forward  scattering  from  these  test  objects.  The 
line  intgrals  of  f(r)  can  then  be  obtained  from 
the  asymptotic  limit  as  k  ♦  •  of  (1.16);  namely, 

lia-rr-v .(r)  -  J  f(r)  <U  (*.2) 

|t*a»  *  L 

where  L  is  the  straight  line  path  from  the  source 
to  our  detector  position.  We  use  the  Thiele  pro¬ 
cedure  described  in  (1.8)  to  (1.14)  to  evaluate 
this  asymptotic  limit  from  an  odd  number  of  fre¬ 
quencies  (usually  9  tqui-epaced  frequencies  from  1 
to  4  Mix). 


We  have  tested  the  application  of  these  meth¬ 
ods  with  s  series  of  computer  simulations  on  the 
example  illustrated  in  Pig.  1.  This  teat  case 
consisted  of  four  objects  with  a  distribution 


given  by 


f(r)  -  e 


These  straight  line  integrals  of  f(r)  were 
eslcuieted  for  e  number  (typically  55)  of  detector 
points  spaced  .3125  mm  apart  along  a  line  oriented 
perpendicular  to  the  direction  of  incident  plane 
wave  end  located  8  cm  from  the  x-y  origin.  Then, 
by  taking  date  fros#  multiple  views  (generally  30 
views  spaced  6*  apart)  and  uaing  ART,  a  wall -known 
X-ray,  computerised  tomography  algorithm,  we  were 


(4.1) 


able  to  reconstruct  the  distribution  of  f(r).  As 
we  have  shown  in  our  previous  paper  (1],  quite 
good  reconstructions  of  f(r)  can  be  obtained  pro* 
vided  we  have  very  little  noise  in  our  date.  An 
illustration  of  a  reconstruction  of  f(r)  from 
reference  1 1)  using  random  noise  with  a  .1  percent 
standard  deviation  is  given  in  Fig.  2.  This  re 


Fig.  2.  Reconstruction  of  the  test  objects  from 
10  sets  of  9  frequencies  (unifomly 
spaced  fro®  1*4  KHz)  with  .1  percent 
random  noise  using  the *1  -norm  after 
applying  the  Thiele  algorithm.  Reprinted 
with  permission  from  reference  1  (Sten* 
ger ,  F . ,  et  al . ,  Proc.  ofConf.  on  Wave 
Motion  held  at  the  University  of  Toronto^ 
July.  1983,  Worth  Holland). 

construction  was  obtained  by  using  10  distinct 
sets  of  9  frequencies  from  1*4  KHz  and  applying 
the  1  *norm  averaging  to  the  results  of  the  Thiele 
algorithm.  However,  the  procedure  of  applying 
Che  2^  averaging  after  the  Thiele  procedure  soon 
fails1 to  produce  acceptable  results  as  we  intro¬ 
duce  increasingly  larger  amounts  of  noise,  since 
the  Thiele  algorithm  magnifies  the  noise.  For 
example.  Fig.  3  gives  s  reconstruct  ion  on  data 
with  1.0  percent  random  noise  using  this  method 
with  54  averages  of  each  Thiele  result.  The  ob* 
jects  are  just  barely  recognizable  end  the  peek 
values  are  low  in  magnitude  by  more  than  a  factor 
of  2x. 

Our  new  procedures  described  in  sections  2 
and  3  can  readily  handle  noiee  of  1  or  2  percent. 
Figure  4  shows  a  reconstruction  of  our  test  case 
from  data  with  1  percent  random  noiae  using  our 
new  procedure  with  54  averages.  Although  all  of 
the  object#  are  quite  well  resolved  in  Fig.  4(  we 
diacovered  that  by  using  a  constraint  condition  in 
our  reconstruction  algorithm  (based  upon  the  fact 
that  f(r)  was  known  to  be  of  one  sign),  me  were 
able  to  even  further  improve  the  resolution,  at 
shown  in  Fig.  5  (notice  the  narrower  line  profiles 
and  how  much  closer  the  troughs  are  to  aero). 
Figure  6  shows  how  much  more  robust  our  new  proce- 


Fig.  3.  Reconstruction  of  the  test  objects  from 
data  with  1  percent  aoiae  using  the  2  * 
norm  procedure  after  applying  the  Thiele 
algorithm,  as  in  the  previous  figure  (but 
here  the  averaging  was  done  with  54  re¬ 
peats  of  noisy  dats  using  9  equally 
spaced  frequencies). 


Fig.  4.  Reconstruction  of  the  test  objects  from 
data  with  1  percent  noiae  ueing  our  new 
2  averaging  and  the  adaptive  feature 
described  in  sections  2  and  3.  As  in 
Fig.  3,  the  ft  averaging  was  obtained  on 
54  repeats  on1  9  equally  spaced  frequen¬ 
cies  from  1*4  MHz. 

dure  it  as  this  reconstruction  was  obtained  from 
data  with  1  percent  noiae  using  only  10  values  in 
our  ft  averaging.  Using  our  previous  method,  we 
were  not  able  to  obtain  recognisable  images  from 
only  10  2^  averages  of  data  with  l  percent  noiae. 

We  also  examined  two  other  interesting  ques¬ 
tions.  The  first  question  concerned  the  applies* 
bility  of  maing  only  the  (2wfc(r)/ik)  term  in 


Fig.  5.  Reconstruction  of  the  test  objects  from 
dsts  with  1  percent  noise  obtained  in  the 
same  manner  as  Fig.  4.,  but  using  a  con¬ 
straint  condition  in  ART  based  upon  know¬ 
ing  that  the  distribution  was  of  one 
sign. 


The  final  question  concerned  the  effect  of 
our  stopping  criteria.  Our  tests  showed  that 
while  the  stopping  procedure  described  in  section 
3  was  quite  helpful,  it  was  not  essential  for  ob¬ 
taining  a  useable  image.  For  example,  Fig.  7 


Fig.  6.  Reconstruction  of  the  test  objects  from 
data  with  1  percent  noise  using  our  new 
procedure  (as  in  Fig.  5),  but  doing  the 
£  averaging  over  only  10  sets  of  data. 
(For  this  reconstruction,  a  constraint 
condition  was  also  incorporated  into  the 
ART  algorithm  triiich  took  advantage  of  the 
knowledge  that  f(r)  was  of  one  sign.) 

(1.16)  at  a  single  wave  number  k  as  a  first  order 
approximation  to  the  asymptotic  limit.  He  learned 
that  for  the  test  object  of  Fig.  1  (FVHM  *  2  n), 
this  first  term  gives  a  completely  unacceptable 
approximation  even  without  any  noise.  He  could 
oot  even  distinguish  the  locations  of  the  objects 
if  we  used  only  this  term. 


Fig.  7.  Reconstruction  of  the  test  objects  from 
data  with  1  percent  noise.  As  for  Fig. 
3,  we  used  the  new  £  averaging  procedure 
on  54  sets  of  noisy  oat  a  before  applying 
the  Thiele  algorithm,  but  here  the  stop¬ 
ping  criteria  waa  omitted  in  the  Thiele 
procedure.  Degradation  from  the  results 
in  Fig.  4  where  the  stopping  criteria  was 
used  is  apparent . 

shows  s  reconstruction  under  the  same  basic  situa¬ 
tion  as  Fig.  4  (i.e.,  9  frequencies  using  £  aver¬ 
aging  54  times  before  applying  the  Thiele  proce¬ 
dure),  but  without  the  stopping  procedure.  While 
the  resolution  in  7  is  much  less  than  Fig.  4,  one 
can  still  clearly  distinguish  all  four  objects. 

5.  Conclusions 

In  a  previous  paper,  we  described  the  basic 
ides  of  transforming  the  ultrasound  inverse  scat¬ 
tering  problem  (with  its  attendant  need  to  treat 
diffraction  effects  in  the  measured  scattering 
data)  to  an  equivalent  X-ray  CT-like  problem  (with 
no  diffraction  effects  in  the  new  transformed 
data)  by  the  uae  of  frequency  extrapolation  [3]. 
The  basic  ides  of  our  method  is  to  remove  diffrac¬ 
tion  effects  by  generating  a  transformed  projec¬ 
tion  of  the  object  at  infinite  acoustic  frequen¬ 
cies  (3).  Our  early  attempts  to  epply  the  method 
met  with  limited  success  because  of  extreme  sensi¬ 
tivity  to  noise  in  the  date.  Applying  £  averag¬ 
ing  to  many  sets  of  frequency  extrapolated  projec¬ 
tions  helped  improve  noise  tolerance,  but  probably 
oot  enough  for  practical  application#  with  real 
data  [1].  Our  latest  algorithm  using  £  averaging 
(median  filtering)  on  the  date  before1  frequency 
extrapolation  combined  with  an  adaptive  frequency 
extrepolet ion  method  provided  even  further  im¬ 
provements.  Images  have  been  made  from  simulated 


data  with  l  percent  noise  with  as  few  »i  10  sets 
of  data  combined  by  t  averaging.  Such  data  could 
be  obtained  with  existing  tomographic  scsnners  at 
frequencies  below  4  MH*.  Thus,  out  newest  algo¬ 
rithm  should  now  be  tested  with  reel  scattering 
date  to  study  its  usefulness  in  ultrasound  dif- 
fraction  tomographic  imaging.  We  also  have  ahown 
in  this  paper  how  to  reconstruct  two  material 
properties  with  differing  dependencies  on  k.  This 
suggests  thst  our  method  might  allow  reconstruc¬ 
tion  of  density  in  addition  to  speed  of  sound; 
however,  this  has  not  been  verified. 
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In  our  previous  work  [1],  a  procedure  was  developed,  based  on  algorithms 
for  rational  function  frequency  extrapolation  and  Hj-norm  averaging  to 
circumvent  both  the  effects  of  noise,  and  errors  due  to  refracted  curved  paths, 
for  solving  the  inverse  scattering  problem  in  ultrasonic  imaging.  In  this 
paper,  we  reverse  the  order  of  the  first  two  algorithms  with  respect  to  [1]; 
that  is,  we  first  apply  the  2,,  averaging  procedure  and  then  carry  out  the 
rational  function  procedure,  extrapolating  to  infinite  frequency.  An  adaptive 
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method  for  choosing  the  best  rational  function  expansion  is  also  employed. 

We  describe  the  underlying  ideas  of  the  procedure,  and  we  illustrate  examples 
of  reconstructions  in  the  presence  and  absence  of  Gaussian  noise.  These 
results  are  then  compared  with  results  of  previous  experiments.  The  resulting 
procedure  works  with  a  much  smaller  signal-to-noise  ratio.  An  extension  is 
suggested  to  image  speed  of  sound  and  density. 
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